%%  ESS 590D: HEAT AND MASS FLOW MODELLING  SPRING 2006
%%
%%  THIS IS A WAVE-EQUATION SOLVER FOR SH WAVES
%%
%%        d^2(rho phi)/dt^2 = div( Gamma grad phi)
%%
%%  phi   = displacement normal to propagation plane
%%  Gamma = elastic shear modulus
%%  rho   = density
%%
%%  This code is adapted from transient diffusion code
%%  The only differences are in the transient term.
%%  Derivation closely follows solution of momentum equations 
%%  to derive Navier Stokes equations.
%%  In elasticity, displacement is the variable that transmits 
%%  momentum into a volume through shear modulus, rather than
%%  velocity transitting momentum through viscosity.
%%  In order to have displacement as the variable inteh transient
%%  term, we have to write the velocity in
%%
%%      [ (rho v)_P - (rho v)_P^0 ] dx dz
%%
%%  as the _rate of change_ of displacement u, giving
%%
%%         [ (rho u_dot)_P - (rho u_dot)_P^0 ] dx dz
%%      =  [ rho( (u_P - u_P^0)/dt - (u_P^0 - u_P^-1)/dt ] dx dz
%%
%%  where u_P^0 and u_P^-1 are the (known) solution at node P at
%%  previous and next previous times.
%%  After dividing through both sides of discretization equation by dt,
%%  we still have a coefficient a_P_0, except that instead of being
%%    a_P^0 = rho dx dz/dt
%%  it is now
%%   a_P^0 = rho dx dz/dt^2
%%  The coefficients a_E, a_W etc are unchanged from the
%%  diffusion problem. 
%%  The right-hand-side term b of known constants,
%%  rather than having [ a_P^0 * u_P^0 ] has instead 
%%      [ 2 * a_P^0 * u_P^0 - a_P^0 * u_P^-1 ]
%%
%%  Caution:  Although this scheme follows all Patankar's rules,
%%            results can be inaccurate if:
%%            - mesh is too coarse to capture shape of wave.
%%            - time step is too large to ensure that temporal curvature
%%              (centered at u_P^0) can balance spatial curvature
%%              (centered on u_P).   
%%              Large temporal mismatch in centers of curvature causes 
%%              dissipation, and waves die out.
%%
%%
%%  Boundary Conditions:
%%   All boundaries are on finite-volume interfaces
%%   (Patankar Practice B, page 69)
%%   However, there is an additional "node" X_B located on each
%    boundary-volume face.
%%   Its value and coefficients are determined by phi value or
%%   gradient boundary conditions there.
%%   Adjacent interior point is X_I
%%   (X_I - X_B) = dX_B
%%
%%    phi_B = phi_0
%%        a_B = 1,  all other a_J = 0
%%        b   = phi_0
%%
%%    d phi/dx = dphi_dx_0
%%        use gradient between X_B and X_I (interior point)
%%        d phi/dx|B = Gamma_I (phi_I - phi_B)/dx_B
%%    or
%%        phi_B = dphi_dx_0 * dx_B + phi_I
%%
%%
%%  vector bc_e  - specified phi gradient on east boundary
%%  vector bc_w  - specified phi gradient on west boundary
%%  vector bc_d  - specified phi gradient on lower boundary
%%  vector bc_u  - specified phi gradient on upper boundary
%%
%%  There must be at least 1 place where a value of phi is set.
%%  This is achieved by setting a phi "source" in the upper boundary.
%%  In the region "source_spot", phi BC is a single cycle of a sine wave
%%  with amplitude "source_height" and period "source_period" in time.
%%  By making "source_spot" longer or shorter it is possible to
%%  simulate everything from plane waves to "point" sources.
%%
%%  Tis code can run 4 different examples, based on the user-set
%%  value of "fig_num"
%%   fig_num = 1   plane wave into homogeneous domain
%%   fig_num = 2   circular (point source) wave in homogeneous domain
%%   fig_num = 3   plane wave into domain with an internal interface
%%   fig_num = 4   plane wave into domain with a fluid inclusion
%%
%%
%%
%%                                   Ed Waddington      May 24 2006
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
   global source_spot  source_period  source_height 

   global dt_plot  fig_num


%
%  seismic source parameters
%    source_spot   - node index on upper boundary for line source
%                    (will be set after grid has been set)
%    source_period - period of single sine-wave cycle
%    source_height - amplitude of displacement source
% -----------------------------------------------------
    source_period = 1;
    source_height = 1;

%
%   fig_num determines nature of model
%   fig_num = 1   plane wave into homogeneous domain
%   fig_num = 2   circular (point source) wave in homogeneous domain
%   fig_num = 3   plane wave into domain with an internal interface
%   fig_num = 4   plane wave into domain with a fluid inclusion
%  ---------------------------------------------------------------
    fig_num = 2;
%

%%%%%%%%%%%%%%%%%%%%%%%%%%
%  set up finite volumes
%%%%%%%%%%%%%%%%%%%%%%%%%%
%   This code sets up uniform-sized volumes (could be generalized)
%   There are points on all boundaries, which are also
%   volume interfaces  (Patankar Practice B, page 69)
%
   Lx    = 2;               %  x-extent of domain
   nx_fv = 50;                %  number of finite volumes in x
   nx_P  = nx_fv + 2;        %  number of nodes in x
%
%
   Lz    = 2;               %  z-extent of domain
   nz_fv = 60;                %  number of finite volumes in z
   nz_P  = nz_fv + 2;        %  number of nodes in z
%
%
%  set up matrices X_edges, Z_edges - defining (x,z) coordinates 
%  of finite-volume edges
%
%    set x-widths of volumes 
%    this bit of code sets up uniform edge spacing
%   ----------------------------------------------
     dx = Lx/nx_fv;              %  z-widths of volumes
     x_edges_vec = dx * (0:nx_fv);
%
  %%   X_edges = repmat(x_edges_vec, nz_fv+1, 1);
%
     dz = Lz/nz_fv;               %  z-widths of volumes
     z_edges_vec = dz * (0:nz_fv)';
%
     Z_edges = repmat(z_edges_vec, 1, nx_fv+1);

%
% set nodes at finite-volume centers and on all boundaries
     x_P_vec = [ x_edges_vec(1) ...
                (x_edges_vec(1:end-1) + diff(x_edges_vec)/2) ...
                 x_edges_vec(end) ]; 
%
     X_P = repmat(x_P_vec, nz_P, 1);
%
%
     z_P_vec = [ z_edges_vec(1) ...
                (z_edges_vec(1:end-1) + diff(z_edges_vec)/2)' ...
                 z_edges_vec(end) ]'; 
%
     Z_P = repmat(z_P_vec, 1, nx_P);
%
%
%
%%  Boundary Conditions 
%%%%%%%%%%%%%%%%%%%%%%%%%
%%
%  Boundary conditions along each boundary are carried in an array
%  labelled bc_x.  For each boundary point, 
%        first entry is numerical value of BC,
%        second entry indicates type of BC
%               1 for specified value
%               2 for specified gradient
%  ---------------------------------------------------
%
%  phi gradient at east boundary
     bc_e_0 = 0;
     bc_type   = 2;
     bc_e = repmat( [ bc_e_0  bc_type ], nz_P, 1 );
%
%%  phi gradient at west boundary
     bc_w_0 = 0;
     bc_type   = 2;
     bc_w = repmat( [ bc_w_0  bc_type ], nz_P, 1 );
%
%
%%  phi gradient on upper boundary
     bc_u_0 = 0;
     bc_type = 2;
     bc_u   = repmat( [ bc_u_0  bc_type ]', 1, nx_P );
%
%  add displacement source on upper boundary at source_spot
  if(fig_num == 1)
 %  plane wave in uniform medium
      source_spot = 1:nx_P;
%
  elseif(fig_num == 2)
 %  "point" source  (3 volumes)
 %%  source_spot = fix(nx_P/2) + 1;
     source_spot = [ (fix(nx_P/2) - 1):(fix(nx_P/2) + 1) ];
%
  elseif(fig_num == 3)
 %  plane wave hits sloping interface
      source_spot = 1:nx_P;
%
  elseif(fig_num == 4)
 %  plane wave hits fluid inclusion
      source_spot = 1:nx_P;
  end
%
     bc_u(2, source_spot) = 1;
%
%
%
%%  phi gradient on lower boundary
     bc_d_0 = 0;
     bc_type = 2;
     bc_d   = repmat( [ bc_d_0  bc_type ]', 1, nx_P );
%

%
%
%  set time-step control parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      dt = 0.002;
      nt = 2000;
%
      dt_plot = 0.1;


%  set up initial condition
%  ------------------------
      phi_0 = zeros( nz_P, nx_P );
%
%  eliminate values on corners where there are no nodes
      bc_e(1,1)    = NaN;
      bc_e(1, end) = NaN;
      bc_w(1,1)    = NaN;
      bc_w(1, end) = NaN;
      bc_u(1,1)    = NaN;
      bc_u(1, end) = NaN;
      bc_d(1,1)    = NaN;
      bc_d(1, end) = NaN;

      phi_0(1,1)     = NaN;
      phi_0(1,end)   = NaN;
      phi_0(end,1)   = NaN;
      phi_0(end,end) = NaN;
%
%%  material parameters
%%%%%%%%%%%%%%%%%%%%%%%
% set up thermal-property matrix - uniform conductivity
     Gamma_0 = 1;                     %   W m^-1 deg^-1    
     Gamma_P = repmat( Gamma_0, nz_P, nx_P );
%
%
%   set up interfaces
  if(fig_num == 1)
 %  plane wave in uniform medium
      %  Gamma_P is set to go;
%
  elseif(fig_num == 2)
 %  "point" source  (3 volumes)
 %%  source_spot = fix(nx_P/2) + 1;
       %  Gamma_P is set to go;
%
  elseif(fig_num == 3)
 %  plane wave hits sloping interface 
 %  slower wave speed below interface
 %   Z = z_0 * (1 - X) + z_1 * X
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
       Gamma_1 = 0.2;
%
       z_range = z_P_vec(end) - z_P_vec(1);
       z_0 = 0.9;
       z_1 = 0.2;
%
       x_range = x_P_vec(end) - x_P_vec(1);
       x_0 = 0;
       x_1 = 1;

       for ix = 1:nx_P
          X_bdry = x_P_vec./ x_range;
          Z_bdry = (z_0 * (1 - X_bdry) + z_1 * X_bdry)';
          index = find( Z_P(:,ix)/z_range >= Z_bdry(ix) ); 
          Gamma_P( index, ix) = Gamma_1;
       end   %  for ix = 1:nx_P ...
%
  elseif(fig_num == 4)
 %  plane wave hits fluid inclusion at some internal points
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   shear modulus of near-fluid inclusion
    Gamma_1 = 0.0001;
    
    x_range = x_P_vec(end) - x_P_vec(1);
    z_range = z_P_vec(end) - z_P_vec(1);

  % center and radius of inclusion
    x_center = x_range/2;
    z_center = z_range/2;
    radius   = x_range/6;
    
  % non-dimensional distance between each node and inclusion center
    R_2 = ( (X_P -x_center)./radius ).^2 + ...
          ( (Z_P -z_center)./radius ).^2;
  %
  %  put Gamma_1 into volumes inside circle
    index = find( R_2 <= 1);
    Gamma_P(index) = Gamma_1;
%  
  end    %  if( fig_num == ...



%
% set up density matrix
     rho_0 = 1;                       %   kg m^-3
     rho = repmat( rho_0, nz_P, nx_P );
%

%
%%  source term   
%%%%%%%%%%%%%%%%%%%%%
%   S = S_C + S_P * phi_P 
%   ----------------------
     S_C_0  = 0;                   %  W m^-3 s^-1
     S_P_0  = 0;                   %  W m^-3 s^-1 deg^-1
%
     S_C  = repmat( S_C_0, nz_P, nx_P );  %  constant source term
     S_P  = repmat( S_P_0, nz_P, nx_P );  %  phi_P-dependent source
%
%

%
% This is it! solve for phi using Finite Volumes
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    phi = fvm_soln_SH_waves_2D( x_edges_vec, z_edges_vec, ...
                  x_P_vec, z_P_vec, nt, dt, Gamma_P, rho, S_C, S_P, ...
                  bc_u, bc_d, bc_e, bc_w, phi_0 );
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
%
%

